Metadata¶

Preprocess cranberry bogs to remove unwanted artifacts and aggregate by proximity

Author: Adrian Wiegman

Date Modified: 10/31/2023

Creates a layer of cranberry bog polygons from raw input data

Main Program¶

In [25]:
# iphython options
# delete variables in workspace
%reset -f
#places plots inline
%matplotlib inline
#automatically reloads modules if they are changed
%load_ext autoreload 
%autoreload 2
# this codeblock sets up the environment from jupyter notebooks
setup_notebook = "C:/Users/Adrian.Wiegman/Documents/GitHub/Wiegman_USDA_ARS/Cran_Q_C/2_gis/scripts/_Setup.ipynb"
%run $setup_notebook # magic command to run the notebook
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
***
loading python modules...

  `module_list` contains names of all loaded modules

...module loading complete

***
loading user defined functions...

type `fn_`+TAB to for autocomplete suggestions

 the object `def_list` contains user defined function names:
   fn_get_info
   fn_arcgis_table_to_df
   fn_arcgis_table_to_np_to_pd_df
   fn_run_script_w_propy_bat
Signature: fn_get_info(name='fn_get_info')
Docstring: returns the source information about a given function name
Source:   
def fn_classify_wetlands(x):
    """
    Consolidates DEP wetlands into fewer categories
    """
    if "MARSH" in x:
        return("MARSH")
    elif "WOOD" in x:
        return("FORESTED SWAMP")
    elif "SHRUB" in x:
        return("SHRUB SWAMP")
    elif len(x) == 0:
        return("NON WETLAND")
    else:
        return(x)
File:      Dynamically generated function. No source code available.
Type:      function
   fn_try_mkdir
   fn_hello
   fn_recursive_glob_search
   fn_regex_search_replace
   fn_regex_search_0
   fn_arcpy_table_to_excel
   fn_agg_sum_df_on_group
   fn_add_prefix_suffix_to_selected_cols
   fn_calc_pct_cover_within_groups
   fn_buildWhereClauseFromList
   fn_FA_to_Q
   fn_alter_field_double
   fn_return_float
   fn_classify_wetlands

 use ??{insert fn name} to inspect
 for example running `??fn_get_info` returns:

 you can also use `fn_get_info(name={insert fn name})` to do the same thing as `??{insert fn name}`

...function loading complete

***
managing project directories:
 working directory (`wdr`): C:\Workspace\Geodata\Cran_Q_C\
 temp dir (`tdr.name`): C:\Users\ADRIAN~1.WIE\AppData\Local\Temp\ArcGISProTemp29724\tmpac5oidy7
 input file directory (`idr`): C:\Workspace\Geodata\Cran_Q_C\inputs
 output dir (`odr`): C:\Workspace\Geodata\Cran_Q_C\outputs




***
setting up arcpy environment...

 default geodatabase path (`gdb_default_path`): C:\Workspace\Geodata\Cran_Q_C\Default.gdb
 active map (`mx`) name: Map2_prep_gw
 output coordinate system: NAD_1983_UTM_Zone_19N

... env setup complete


+++SETUP COMPLETE+++


In [26]:
#wdr = r"C:\Workspace\Geodata\Cran_Q_C"
print(wdr)
C:\Workspace\Geodata\Cran_Q_C\
In [27]:
# Select file geodatabase to store results\n",
#gdb = "Cran_Q_C.gdb"
#ap.env.workspace = os.path.join(wdr,gdb)
In [28]:
# copy raw bog data to the default geodatabase in working directory
in_feat = r"C:\Workspace\Geodata\Cran_Q_C\inputs\WMAbogsDRAFT2013\WMAbogsDRAFT2013.shp"
bogs_raw = "bogs_raw"
arcpy.management.CopyFeatures(
    in_features=in_feat,
    out_feature_class=bogs_raw)
Out[28]:

Messages

Start Time: Wednesday, December 4, 2024 10:04:57 AM
Succeeded at Wednesday, December 4, 2024 10:04:58 AM (Elapsed Time: 1.14 seconds)

process cranberry bogs to split features.¶

Example of unwanted artifact of cranberry bog layer.

Polygons that have the same owner are joined together by thin corridors. This can be problematic if the bog overlaps with a stream reach (denoted by the green point). These corridors sections must be removed. The corridors can be removed by putting a negative buffer on the the polygons. Additionally it may be helpful to split polygons that are far apart. And to combine polygons that are close together.

image-2.png

In [29]:
# Negative buffer around perimeter to remove small corridors (this removes links bewteen bogs with the same permit that are not hydrologically connected)
buff_dist = 3
bogs_neg = "bogs_buff_neg{}m".format(buff_dist)
arcpy.analysis.PairwiseBuffer(
    in_features=bogs_raw,
    out_feature_class=bogs_neg,
    buffer_distance_or_field="-{} Meters".format(buff_dist),
    dissolve_option="NONE",
    dissolve_field=None,
    method="PLANAR",
    max_deviation="0 Meters"
)
---------------------------------------------------------------------------
ExecuteError                              Traceback (most recent call last)
Cell In[29], line 4
      2 buff_dist = 3
      3 bogs_neg = "bogs_buff_neg{}m".format(buff_dist)
----> 4 arcpy.analysis.PairwiseBuffer(
      5     in_features=bogs_raw,
      6     out_feature_class=bogs_neg,
      7     buffer_distance_or_field="-{} Meters".format(buff_dist),
      8     dissolve_option="NONE",
      9     dissolve_field=None,
     10     method="PLANAR",
     11     max_deviation="0 Meters"
     12 )

File C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py:1298, in PairwiseBuffer(in_features, out_feature_class, buffer_distance_or_field, dissolve_option, dissolve_field, method, max_deviation)
   1296     return retval
   1297 except Exception as e:
-> 1298     raise e

File C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py:1281, in PairwiseBuffer(in_features, out_feature_class, buffer_distance_or_field, dissolve_option, dissolve_field, method, max_deviation)
   1277 from arcpy.arcobjects.arcobjectconversion import convertArcObjectToPythonObject
   1279 try:
   1280     retval = convertArcObjectToPythonObject(
-> 1281         gp.PairwiseBuffer_analysis(
   1282             *gp_fixargs(
   1283                 (
   1284                     in_features,
   1285                     out_feature_class,
   1286                     buffer_distance_or_field,
   1287                     dissolve_option,
   1288                     dissolve_field,
   1289                     method,
   1290                     max_deviation,
   1291                 ),
   1292                 True,
   1293             )
   1294         )
   1295     )
   1296     return retval
   1297 except Exception as e:

File C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py:532, in Geoprocessor.__getattr__.<locals>.<lambda>(*args)
    530 val = getattr(self._gp, attr)
    531 if callable(val):
--> 532     return lambda *args: val(*gp_fixargs(args, True))
    533 else:
    534     return convertArcObjectToPythonObject(val)

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Feature Layer.
WARNING 000725: Output Feature Class: Dataset C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_buff_neg3m already exists.
Failed to execute (PairwiseBuffer).
In [ ]:
# positive buffer to return to original area
bogs_neg_pos = bogs_neg+"_pos{}m".format(buff_dist)
arcpy.analysis.PairwiseBuffer(
    in_features=bogs_neg,
    out_feature_class=bogs_neg_pos,
    buffer_distance_or_field="{} Meters".format(buff_dist),
    dissolve_option="NONE",
    dissolve_field=None,
    method="PLANAR",
    max_deviation="0 Meters"
)

Screenshot%202023-10-23%20134356.png

In [ ]:
# split all features containing multipart polygons (unconnected) into single part features (one feature per polygon). 
bogs_split = bogs_neg_pos+"_split"
arcpy.management.MultipartToSinglepart(bogs_neg_pos, bogs_split)
In [ ]:
# aggregate polygons of bogs within 30 meters and sharing the same FID. 
# note that this modifies the shape of polygons by creating new corridors that link between adjacent polygons
# this increases total area to the of all polygons in the layer (which can lead to overestimates of potentially restorable bog area). 
agg_dist = 30 # meters
bogs_agg = "bogs_agg"
arcpy.cartography.AggregatePolygons(
    in_features=bogs_split,
    out_feature_class=bogs_agg,
    aggregation_distance="{} Meters".format(agg_dist),
    minimum_area=None,
    minimum_hole_size="0 SquareMeters",
    orthogonality_option="NON_ORTHOGONAL",
    barrier_features=None,
    out_table="{}_Tbl".format(bogs_agg),
    aggregate_field="ORIG_FID"
)
In [ ]:
# use identity to get FID of aggregated bogs 
bogs_ident = "bogs_ident"
arcpy.analysis.Identity(
    in_features=bogs_split,
    identity_features=bogs_agg,
    out_feature_class=bogs_ident,
    join_attributes="ONLY_FID",
    cluster_tolerance=None,
    relationship="NO_RELATIONSHIPS"
)
In [ ]:
# dissolve all features with that have the same value for FID of bog_agg
# this groups feature within 30m but without adding area 
bogs_ident_diss = "bogs_ident_dissolved"
arcpy.management.Dissolve(
    in_features=bogs_ident,
    out_feature_class=bogs_ident_diss,
    dissolve_field="FID_{}".format(bogs_agg),
    statistics_fields=None,
    multi_part="MULTI_PART",
    unsplit_lines="DISSOLVE_LINES",
    concatenation_separator=""
)
In [ ]:
# use identity to get original raw bog FIDs
bogs_ident_diss = "bogs_ident_dissolved"
bogs_raw = "bogs_raw"
bogs_agg_diss_ident = "bogs_agg_diss_ident"
arcpy.analysis.Identity(
    in_features=bogs_ident_diss,
    identity_features=bogs_raw,
    out_feature_class=bogs_agg_diss_ident,
    #join_attributes="ONLY_FID",
    cluster_tolerance=None,
    relationship="NO_RELATIONSHIPS")
In [ ]:
# Add geometry attributes to table
arcpy.management.CalculateGeometryAttributes(
    in_features="bogs_agg_diss_ident",
    geometry_property="bogsurf_area_m2 AREA;bogsurf_perim_m2 PERIMETER_LENGTH",
    length_unit="METERS",
    area_unit="SQUARE_METERS",
    coordinate_system=None,
    coordinate_format="SAME_AS_INPUT")
In [ ]:
# Remove unwanted artifacts of preprocessing. 


# select and delete any polygons that do not overlap with original bogs 
# these areas are noted by FID_bogs_raw = -1
bogs_agg_diss_ident = "bogs_agg_diss_ident"
arcpy.management.SelectLayerByAttribute(
    in_layer_or_view=bogs_agg_diss_ident,
    selection_type="NEW_SELECTION",
    where_clause="FID_bogs_raw = -1",
    invert_where_clause=None)
In [ ]:
arcpy.management.DeleteRows(
    in_rows=bogs_agg_diss_ident)
In [ ]:
bogs_agg_diss_ident = "bogs_agg_diss_ident"
arcpy.management.SelectLayerByAttribute(
    in_layer_or_view=bogs_agg_diss_ident,
    selection_type="CLEAR_SELECTION")

# remove rows where area is <100m
arcpy.management.SelectLayerByAttribute(
    in_layer_or_view=bogs_agg_diss_ident,
    selection_type="NEW_SELECTION",
    where_clause="Shape_Area < 100")
In [ ]:
arcpy.management.DeleteRows(in_rows=bogs_agg_diss_ident)

pink areas show farm areas after splitting but before aggregation red areas with black outline show area after aggregation green points are the points of maximum flow accumulation within each aggregated cluster of beds

image-2.png

In [ ]:
copy = arcpy.management.CopyFeatures(
    in_features=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_raw",
    out_feature_class=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_raw_copy",
    config_keyword="",
    spatial_grid_1=None,
    spatial_grid_2=None,
    spatial_grid_3=None)
arcpy.management.CalculateField(
    in_table=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_raw_copy",
    field="FID_bogs_raw",
    expression="!OBJECTID!",
    expression_type="PYTHON3",
    code_block="",
    field_type="",
    enforce_domains="NO_ENFORCE_DOMAINS")
In [ ]:
copy = arcpy.management.CopyFeatures(
    in_features=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_agg",
    out_feature_class=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_agg_raw",
    config_keyword="",
    spatial_grid_1=None,
    spatial_grid_2=None,
    spatial_grid_3=None)
In [ ]:
arcpy.management.JoinField(
    in_data="bogs_agg_raw",
    in_field="ORIG_FID",
    join_table=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_raw",
    join_field="OBJECTID",
    fields=None,
    fm_option="NOT_USE_FM",
    field_mapping=None,
    index_join_fields="NO_INDEXES")
In [ ]:
arcpy.analysis.SummarizeWithin(
    in_polygons=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_agg_raw",
    in_sum_features="bogs_agg_diss_ident",
    out_feature_class=r"C:\Workspace\Geodata\Cran_Q_C\Default.gdb\bogs_agg_raw_SummarizeWithin",
    keep_all_polygons="KEEP_ALL",
    sum_fields="bogsurf_area_m2 Sum;bogsurf_perim_m2 Sum",
    sum_shape="NO_SHAPE_SUM",
    shape_unit="SQUAREKILOMETERS",
    group_field=None,
    add_min_maj="NO_MIN_MAJ",
    add_group_percent="NO_PERCENT",
    out_group_table=None)